#include "stdafx.h"
#include "CParticleSystem.h"
#include "graphics/OpenGLWnd.h"
#include "graphics/Texture.h"
#include <vector>
#include <string>

#define _USE_MATH_DEFINES
#include <math.h>


#ifdef _DEBUG
#define new DEBUG_NEW
#endif

void drawCorner(float *cor, float *L, float *R, float ht, int dir);
typedef std::vector<CParticleSystem::CParticle>::iterator pIterator;
float randFloat(float low=0.f, float high=1.f);


//----------------------------------------Constructors------------------------------------------
//  CParticleSystem
//  CParticleSystem::CParticle
//  initParams
//----------------------------------------------------------------------------------------------

CParticleSystem::CParticleSystem(int x, int y, int z, int pcfg)
{
	pConfig = pcfg;

	wood.LoadFile(L"graphics/wood_tex.bmp");
	cy = 0.;
	maxX = x; //maxX = 12;
	maxY = y; //maxY = 12;
	maxZ = z; //maxZ = 12;
	debug = false;
	add = false;
		
	//For CParticle::updateColor(), set a maximum expected velocity. Might need adjustment.
	Vector3df maxV(2.f, 2.f, 2.f);
	maxVecLen = sqrt(maxV.dot(maxV));

	if (pConfig == 2)
	{
		float vars[] = {1.f, 2.f, 6.f, 11.f};
		initGrid(vars, .3f);
	}
	else if (pConfig == 1)
	{
		float vars[] = {3.f, 9.f, 9.f, 11.f};
		initGrid(vars, .6f);
	}
	else if (pConfig == 3)
	{
		float vars[] = {1.f, 2.f, 6.f, 11.f};
		initGrid(vars, .3f);
	}

	initParams();
}

void CParticleSystem::initGrid(float *pars, float del)
{
	//The x,y,z might be in the wrong order... caution!
	index = 0;//initialize the starting index number for new particles

	//Initializing the grid. If we want a block of particles at the top of the grid,
	//  keep using this method, but check the comments when adding a new cell to p2
	if (pConfig == 1 || pConfig == 2) {
		float delta = del;
		for(float curX = pars[0]; curX < pars[1]; curX += delta)
		{
			for(float curY = pars[2]; curY < pars[3]; curY += delta)
			{
				for(float curZ = pars[0]; curZ < pars[1]; curZ += delta)
				{
					particles.push_back(CParticle(index, Vector3df(curX, curY, curZ)));
					index += 1;
				}
			}
		}
	}

	else if (pConfig == 3) 
	{
		//add bottom fluid particles
		float gap = 0.5f;
		for (float botX = 0.f; botX < maxX; botX += gap) 
		{
			for (float botY = 9.f; botY < 12.f; botY += gap)
			//for (float botY = 1.f; botY < 4.f; botY += gap)
			{
				for (float botZ = 0.f; botZ < maxZ; botZ += gap) 
				{
					particles.push_back(CParticle(index, Vector3df(botX, botY, botZ)));
					index += 1;
				}
			}
		}

	}

	// initialize the force vectors
	externalForces.resize(index, Vector3df(0,0,0));

	// initilize the size of grid
	grid.resize(maxX);
	for(int i = 0; i < maxX; i++) 
	{
		grid[i].resize(maxY);
		for(int j = 0; j < maxY; j++)
		{
			grid[i][j].resize(maxZ);
		}
	}

	accelerations.assign(index, Vector3df(0.f, 0.f, 0.f));
}

CParticleSystem::CParticleSystem()
{
	(*this) = CParticleSystem(9,9,9,1);
}

CParticleSystem::CParticleSystem(int p)
{
}

CParticleSystem::~CParticleSystem()
{
}

//Particle Constructors
CParticleSystem::CParticle::CParticle(int i, Vector3df gridPoint)
{
	index = i;

	//pos.setX(gridPoint.x() + randFloat(0.f, 0.99f));
	//pos.setY(gridPoint.y() + randFloat(0.f, 0.99f));
	//pos.setZ(gridPoint.z() + randFloat(0.f, 0.99f));
	pos = gridPoint;

	c[0] = .2f;
	c[1] = .2f;
	c[2] = .8f;
}

CParticleSystem::CParticle::CParticle(CParticle* p)
{
	index = p->index;
	pos = p->pos;
	vel = p->vel;
	c[0] = p->c[0];
	c[0] = p->c[0];
	c[0] = p->c[0];
}

void CParticleSystem::initParams() 
{
	pRadius = .2f;
	leapFrog = true;
	interaction = false;
	externalForce = Vector3df(0.f, 0.f, 0.f);
	
	// ----------------------------------------------------------------------
	// STABLE PARAMETERS (for 3x3x9 box)
	gravity = Vector3df(0.f, -1.f, 0.f);
	MASS = 200.f;		// mass 
	damping = 0.1f;		// damping coefficient
	h = 1.f;			// smoothing length
	k = -20.f;			// gas constant
	rho0 = 500.f;		// rest density
	mew = 50.f;			// fluid viscosity
	deltaT = 0.1f;		// time step
	// ----------------------------------------------------------------------

	if (pConfig == 1)
	{
		gravity = Vector3df(0.f, -G_MAG, 0.f);
		//sigma = 0.6f;		// surface tension coefficient
		//threshold = .2f;	// threshold at which surface tension gets evaluated
		sigma = 0.6f;
		threshold = 1.f;
		deltaT = 0.1f;
	}
	else if (pConfig == 2)
	{
		//gravity = Vector3df(0.f, -1.f, 0.f);
		gravity = Vector3df(0.f, -9.81f, 0.f);
		sigma = 0.6f;		// surface tension coefficient
		threshold = .4f;	// threshold at which surface tension gets evaluated
	}
	else if (pConfig == 3)
	{
		gravity = Vector3df(0.f, -9.81f, 0.f);
		sigma = .6f;		// surface tension coefficient
		threshold = .4f;	// threshold at which surface tension gets evaluated
		deltaT = 0.1f;
	}

	// ----------------------------------------------------------------------
	// STABLE PARAMETERS (for 12x12x12 box):
	//gravity = Vector3df(0.f, -G_MAG, 0.f);
	//MASS = 200.f;			// mass 
	//damping = 0.1f;		// damping coefficient
	//h = 1.f;				// smoothing length
	//k = -20.f;			// gas constant
	//rho0 = 500.f;			// rest density
	//mew = 50.f;			// fluid viscosity
	//sigma = 0.f;			// surface tension coefficient - 0 = no surface tension forces
	//threshold = .2f;		// threshold at which surface tension gets evaluated
	//deltaT = 0.1f;		// time step
	// ----------------------------------------------------------------------


}






//---------------------------------Simulate/Update Functions--------------------------------------
//  Simulate
//  handleCollisions
//  initSimulationStep
//  updateGrid
//  findNeighbors
//  putNeighbors
//------------------------------------------------------------------------------------------------


void CParticleSystem::Simulate()
{
	if (debug)
		int nothing = 0;


	initSimulationStep();
	
	// Step 0: sort particles into grid cells; find neighbors of each particle
	updateGrid();
	findNeighbors();
 
	//// Step 1: compute the density at particle locations 
	Vector3df r;
	for (std::vector<CParticle>::iterator pd = particles.begin(); pd < particles.end(); pd++) {
		CParticle *p = &(*pd);
		if ((int)neighbors[p->index].size() > 0) {
			int pNum = p->index;
			int numNeighbors = (int)neighbors[p->index].size();
		}
		for (std::vector<CParticle>::iterator nd = neighbors[p->index].begin(); nd < neighbors[p->index].end(); nd++) {
			CParticle *n = &(*nd);
			Vector3df r = p->pos - n->pos;
			if (r.length() <= h) {
				float incDensity = MASS * W_poly6(r,h);
				int i = p->index;
				densities[i] += incDensity;
				//densities[pd->index] += MASS * W_poly6(r, h);
			}
			else {
				float nothing = 0.f;
			}
		}
	}

	// Step 2: compute pressure force, viscosity force, and color field for each particle
	float density_p, density_n;
	float pressure_p, pressure_n;
	for (std::vector<CParticle>::iterator pf = particles.begin(); pf < particles.end(); pf++) {
		CParticle *p = &(*pf);
		for (std::vector<CParticle>::iterator nf = neighbors[pf->index].begin(); nf < neighbors[pf->index].end(); nf++) {
			CParticle *n = &(*nf);
			r = p->pos - n->pos;
			if (r.length() <= h) {
				density_p = densities[p->index];
				density_n = densities[n->index];

				pressure_p = k * (density_p - rho0);
				pressure_n = k * (density_n - rho0);

				pressureForces[p->index] +=  MASS * (pressure_p + pressure_n) / (2.f * density_n) * gradientW_spiky(r, h);
				
				viscosityForces[p->index] += mew * MASS * (n->vel - p->vel) / density_n * laplacianW_viscosity(r, h);
				if (viscosityForces[p->index].length() > 0.f) {
					float test = -999.f;
				}

				colorField_gradient[p->index] += MASS / density_n * gradientW_poly6(r, h);
			
				colorField_laplacian[p->index] += MASS / density_n * laplacianW_poly6(r, h);
			}
		}
	}

	// Step 3: move the particles
	Vector3df surfaceTension, totalForce, acceleration;
	for (std::vector<CParticle>::iterator pm = particles.begin(); pm < particles.end(); pm++) {
		CParticle *p = &(*pm);
		// calculate surface tension force
		if (colorField_gradient[p->index].length() >= threshold) {
			surfaceTension = -sigma * colorField_laplacian[p->index] * colorField_gradient[p->index] / colorField_gradient[p->index].length();
		}
		else {
			surfaceTension = Vector3df(0.f, 0.f, 0.f);
		}

		totalForce = surfaceTension + pressureForces[p->index] + viscosityForces[p->index];

		// NOTE: Which one do we use here? The second one looks more realistic, but makes the simulation appear 
		// slower (with the current "stable" parameters we are using.
		//acceleration = totalForce / densities[p->index] * deltaT + gravity + externalForce;
		acceleration = totalForce / densities[p->index] * deltaT + gravity * deltaT + externalForce;

		if (leapFrog) {
			// Leap-Frog integration:
			p->pos += p->vel * deltaT + 0.5f * accelerations[p->index] * pow(deltaT,2);
			p->vel += damping * (0.5f * (accelerations[p->index] + acceleration) * deltaT);
		}
		else {
			// Euler integration:
			p->pos += p->vel * deltaT;
			p->vel += damping * acceleration * deltaT;
		}

		accelerations[p->index] = acceleration;	// Leap-Frog needs accelerations of previous time step
	}

	// Step 4: collision detection
	handleCollisions();
	interaction = false;
}

void CParticleSystem::handleCollisions()
{
	// TO DO: try modeling boundaries as elastic springs with high spring constants
	// like in: http://web.mit.edu/alerer/www/sph.pdf
	// This might help with motion dampening.

	// collision detection
	for(pIterator pi = particles.begin(); pi < particles.end(); pi++)
	{
		CParticle *p = &(*pi);
		float curX = p->pos.x();
		float curY = p->pos.y();
		float curZ = p->pos.z();

		if(curX > maxX || curX < 0)
		{
			p->vel.setX(p->vel.x() * -.9f);
			float diff = curX - (float)maxX;
			if(diff > 0) p->pos.setX((float)maxX - fmod(curX,(float)maxX));
			else p->pos.setX(abs(fmod(curX,(float)maxX)));
		}

		if(curY > maxY || curY < 0)
		{
			p->vel.setY(p->vel.y() * -.75f);
			float diff = curY - (float)maxY;
			if(diff > 0) p->pos.setY((float)maxY - fmod(curY,(float)maxY));
			else p->pos.setY(abs(fmod(curY,(float)maxY)));
		}

		if(curZ > maxZ || curZ < 0)
		{
			p->vel.setZ(p->vel.z() * -.9f);
			float diff = curZ - (float)maxZ;
			if(diff > 0) p->pos.setZ((float)maxZ - fmod(curZ,(float)maxZ));
			else p->pos.setZ(abs(fmod(curZ,(float)maxZ)));
		}
	}
}

void CParticleSystem::initSimulationStep()
{
	if (add)
		addParticles();

	// clear and initialize all
	densities.assign(index,0.f);
	colorField_gradient.assign(index, Vector3df(0.f, 0.f, 0.f));
	colorField_laplacian.assign(index, 0.f);
	pressureForces.assign(index, Vector3df(0.f, 0.f, 0.f));
	viscosityForces.assign(index, Vector3df(0.f, 0.f, 0.f));
	externalForces.assign(index, Vector3df(0.f, 0.f, 0.f));
	accelerations.resize(index);
	neighbors.clear();

	// clear grid
	for(int x = 0; x < maxX; x++)
		for(int y = 0; y < maxY; y++)
			for(int z = 0; z < maxZ; z++)
				grid[x][y][z].pcs.clear();

	if (!interaction)
		externalForce = Vector3df(0.f, 0.f, 0.f);

}

// Sort all particles into proper grid cells using their current positions
void CParticleSystem::updateGrid()
{
	std::vector<CParticle>::iterator pi;
	for(pi = particles.begin(); pi < particles.end(); pi++)
	{
		// current particles position
		Vector3df *pp = &(*pi).pos;

		// NOTE: this "if" statement makes the system not crash, but is NOT the proper way we should be handling these particles
		// that were so far out of the box that collision detection simply projected them to another location out of the box.
		if ((int)pp->x() < 0 || (int)pp->x() >= maxX || (int)pp->y() < 0 || (int)pp->y() >= maxY || (int)pp->z() < 0 || (int)pp->z() >= maxZ) {
			std::string blah = "Houston, we have a problem!";
		}
		else {
			grid[(int)pp->x()][(int)pp->y()][(int)pp->z()].pcs.push_back(&(*pi));
		}
	}
}

// Find potential neighbors of all particles using the grid. Make use of symmetry to improve efficiency.
void CParticleSystem::findNeighbors()
{
	for (int p = 0; p < index; p++) {
		std::vector<CParticle> n; 
		neighbors.push_back(n);
	}
	
	for (int x = 0; x < (int)grid.size(); x++) {
		for (int y = 0; y < (int)grid[0].size(); y++) {
			for (int z = 0; z < (int)grid[0][0].size(); z++) {
				// if grid cell is not empty
				if ((int)grid[x][y][z].pcs.size() != 0) {
					// all particles in own cell are neighbors
					putNeighbors(Vector3df(x, y, z), Vector3df(x, y, z));
					// level y-1
					if (y-1 >= 0) {
						putNeighbors(Vector3df(x, y, z), Vector3df(x, y-1, z));
						if (x-1 >= 0) {
							putNeighbors(Vector3df(x, y, z), Vector3df(x-1, y-1, z));
						}
						if (z-1 >= 0) {
							putNeighbors(Vector3df(x, y, z), Vector3df(x, y-1, z-1));
							if (x-1 >= 0) {
								putNeighbors(Vector3df(x, y, z), Vector3df(x-1, y-1, z-1)); 
							}
							if (x+1 < (int)grid.size()) {
								putNeighbors(Vector3df(x, y, z), Vector3df(x+1, y-1, z-1));
							}
						}
					}
					// level y:
					if (x-1 >= 0) {
						putNeighbors(Vector3df(x, y, z), Vector3df(x-1, y, z));
					}
					if (z-1 >= 0) {
						putNeighbors(Vector3df(x, y, z), Vector3df(x, y, z-1));
						if (x-1 >= 0) {
							putNeighbors(Vector3df(x, y, z), Vector3df(x-1, y, z-1));
						}
						if (x+1 < (int)grid.size()) {
							putNeighbors(Vector3df(x, y, z), Vector3df(x+1, y, z-1));
						}
					}

					// level y+1
					if (y+1 < (int)grid[0].size()) {
						if (x-1 >= 0) {
							putNeighbors(Vector3df(x, y, z), Vector3df(x-1, y+1, z));
						}
						if (z-1 >= 0) {
							putNeighbors(Vector3df(x, y, z), Vector3df(x, y+1, z-1));
							if (x-1 >= 0) {
								putNeighbors(Vector3df(x, y, z), Vector3df(x-1, y+1, z-1));
							}
							if (x+1 < (int)grid.size()) {
								putNeighbors(Vector3df(x, y, z), Vector3df(x+1, y+1, z-1));
							}
						}
					}

				} //if cell not empty
			} //for z
		} //for y
	} //for x
}

// Parameters to and from are the 3D indices of two grid cells  
void CParticleSystem::putNeighbors(Vector3df to, Vector3df from) {

	int tx, ty, tz; int fx, fy, fz;
	tx = (int)to.x(); ty = (int)to.y(); tz = (int)to.z();
	fx = (int)from.x(); fy = (int)from.y(); fz = (int)from.z();

	// if to and from are same grid cell
	if (to == from) {
		for (int f = 0; f < (int)grid[fx][fy][fz].pcs.size(); f++) {
			for (int t = 0; t < (int)grid[tx][ty][tz].pcs.size(); t++) {
				// don't want a particle to be its own neighbor
				if (grid[fx][fy][fz].pcs[f]->index != grid[tx][ty][tz].pcs[t]->index) {
					// do both ways
					CParticle nt = grid[tx][ty][tz].pcs[t];
					CParticle nf = grid[fx][fy][fz].pcs[f];
					if ((nt.pos - nf.pos).length() <= h) {
						neighbors[grid[fx][fy][fz].pcs[f]->index].push_back(nt);
						neighbors[grid[tx][ty][tz].pcs[t]->index].push_back(nf);
					}
				}
			}
		}
	}
	// if to and from are not same grid cell
	else {
		for (int f = 0; f < (int)grid[fx][fy][fz].pcs.size(); f++) {
			for (int t = 0; t < (int)grid[tx][ty][tz].pcs.size(); t++) {
				// do both ways from->to and to->from
				CParticle nt = grid[tx][ty][tz].pcs[t];
				CParticle nf = grid[fx][fy][fz].pcs[f];
				if ((nt.pos - nf.pos).length() <= h) {
					neighbors[grid[fx][fy][fz].pcs[f]->index].push_back(nt);
					neighbors[grid[tx][ty][tz].pcs[t]->index].push_back(nf);
				}
			}
		}
	}
}

void CParticleSystem::addParticles()
{
	if (pConfig == 3)
	{
		float gap = 0.4f;
		for (float curX = 0.f; curX < 1.f; curX += gap)
		{
			for (float curY = 9.f; curY < 12.f; curY += gap)
			{
				for (float curZ = 0.f; curZ < 1.f; curZ += gap) 
				{
					particles.push_back(CParticle(index, Vector3df(curX, curY, curZ)));
					index += 1;
				}
			}
		}
	}
	add = false;
}






//------------------------------------Force Calculations------------------------------------------
//  W_poly6
//  gradientW_poly6
//  laplacianW_poly6
//  gradientW_spiky
//  laplacianW_viscosity
//------------------------------------------------------------------------------------------------

float CParticleSystem::W_poly6(Vector3df r, float h)
{
	float poly6 = (float)(315.f / (64.f *M_PI*pow(h,9))) * pow((pow(h,2) - pow(r.length(),2)),3);
	return poly6;
}

Vector3df CParticleSystem::gradientW_poly6(Vector3df r, float h)
{
	Vector3df gpoly6 = -r * (float)(945.f /(32.f *M_PI*pow(h,9))) * pow((pow(h,2)-pow(r.length(),2)),2);
	return gpoly6;
}

float CParticleSystem::laplacianW_poly6(Vector3df r, float h)
{
	float lpoly6 = (float)(945.f / (8.f *M_PI*pow(h,9))) * (float)(pow(h,2)-pow(r.length(),2)) * (float)(pow(r.length(),2)-3./4.*(pow(h,2)-pow(r.length(),2)));
	return lpoly6;
}

Vector3df CParticleSystem::gradientW_spiky(Vector3df r, float h)
{
	Vector3df gspiky = -r * (float)(45./(M_PI*pow(h,6)*r.length())) * pow(h-r.length(),2);
	return gspiky;
}

float CParticleSystem::laplacianW_viscosity(Vector3df r, float h)
{
	float lviscosity = (float)(45. / (M_PI * pow(h, 5))) * (float)(1. - r.length()/h);
	return lviscosity;
}





//------------------------------------Drawing Functions------------------------------------------
//  Draw
//  drawParticles
//     CParticle::updateColor
//  drawTable
//  drawGlass
//-----------------------------------------------------------------------------------------------

void CParticleSystem::Draw()
{
	glEnable(GL_LIGHTING);
	glEnable(GL_COLOR_MATERIAL);
	//Draw the environment
	drawTable();

	//Draw the particles
	drawParticles();

	//Draw the glass afterwards (for transparency)
	drawGlass();

}


void CParticleSystem::drawParticles()
{
	glEnable(GL_LIGHTING);

	for(pIterator pi = particles.begin(); pi < particles.end(); pi++)
	{
		(*pi).updateColor(maxVecLen);
		drawSphereAt((*pi).pos, (*pi).c, false);
	}


	glDisable(GL_LIGHTING);
}



//Updates particle colors based on velocities
void CParticleSystem::CParticle::updateColor(float maxLen)
{
	//The idea behind this is to convert the length of the current particle's velocity vector
	//  Then, based on an estimation of the greatest velocity we can expect in the system,
	//  creates a hue value with red being the longest velocity and blue the shortest 
	//      (think of a temperature map: red = hot, blue = cold)

	//HSL conversion help from http://130.113.54.154/~monger/hsl-rgb.html
	float curLen = sqrt(vel.dot(vel));
	if (curLen/maxLen > 2.)
	{
		c[0] = .8f;
		c[1] = .2f;
		c[2] = .8f;
		return;
	}

	//Get proportion of current vel / max vel
	//max of 240 (blue), reverses so high = red, low = blue
	float p = 240.f - ((curLen/maxLen)*240.f);
	float H = p/360.f;
	float S = 1.f;
	float L = 0.5f;

	float t2 = L+S - L*S;
	float t1 = 2.f * L - t2;

	float nR, nG, nB;
	float tR = H + 1.f/3.f;
	float tG = H;
	float tB = H - 1.f/3.f;
	

	//Look at the website for more info on these
	if (tR < 0) tR += 1.0f;
	else if (tR > 1) tR -= 1.0f;
	if (tG < 0) tG += 1.0f;
	else if (tG > 1) tG -= 1.0f;
	if (tB < 0) tB += 1.0f;
	else if (tB > 1) tB -= 1.0f;
	
	//Set new red
	if (6.f * tR < 1) nR = t1 + (t2-t1) * 6.f * tR;
	else if (2.f * tR < 1) nR = t2;
	else if (3.f * tR < 2) nR = t1 + (t2-t1) * ((2.f/3.f)-tR)*6.f;
	else nR = t1;

	//Set new green
	if (6.f * tG < 1) nG = t1 + (t2-t1) * 6.f * tG;
	else if (2.f * tG < 1) nG = t2;
	else if (3.f * tG < 2) nG = t1 + (t2-t1) * ((2.f/3.f)-tG)*6.f;
	else nG = t1;

	//Set new blue
	if (6.f * tB < 1) nB = t1 + (t2-t1) * 6.f * tB;
	else if (2.f * tB < 1) nB = t2;
	else if (3.f * tB < 2) nB = t1 + (t2-t1) * ((2.f/3.f)-tB)*6.f;
	else nB = t1;

	c[0] = 100*(nR/255); //Convert to 0-100 range, not 0-255
	c[1] = 100*(nG/255);
	c[2] = 100*(nB/255);
}




//Drawing the simple table on which to place the particle system
//   Center of the surface of the table is at 0,cy,0 for simplicity
//	 cy = member variable
void CParticleSystem::drawTable()
{
	double bx = maxX/2;
	//double by = maxY/2;
	double bz = maxZ/2;
	double wd = 15.;
	double ht = 1.;
	double dp = 15.;

	double ny = cy - .2; // Lower the table!

	GLdouble a[] = {-wd+bx, ny, -dp+bz};
    GLdouble b[] = {wd+bx, ny, -dp+bz};
    GLdouble c[] = {wd+bx, ny, dp+bz};
    GLdouble d[] = {-wd+bx, ny, dp+bz};
    GLdouble e[] = {-wd+bx, ny-ht, -dp+bz};
    GLdouble f[] = {wd+bx, ny-ht, -dp+bz};
    GLdouble g[] = {wd+bx, ny-ht, dp+bz};
    GLdouble h[] = {-wd+bx, ny-ht, dp+bz};
	
    glColor3d(1., 1., 1.);

	glEnable(GL_TEXTURE_2D);
    glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
    glBindTexture(GL_TEXTURE_2D, wood.TexName());
	// Top
    glBegin(GL_QUADS);
		glNormal3d(0, 0, 1);
		glTexCoord2f(0, 0);
		glVertex3dv(d);
		glTexCoord2f(1, 0);
		glVertex3dv(c);
		glTexCoord2f(1, 1);
		glVertex3dv(b);
		glTexCoord2f(0, 1);
		glVertex3dv(a);
    glEnd();

	
    glColor3d(.8, .8, .8);
	// Front
    glBegin(GL_QUADS);
		glNormal3d(0, 0, 1);
		glTexCoord2f(0, 0);
		glVertex3dv(d);
		glTexCoord2f(1, 0);
		glVertex3dv(h);
		glTexCoord2f(1, 1);
		glVertex3dv(g);
		glTexCoord2f(0, 1);
		glVertex3dv(c);
    glEnd();

	glColor3d(.4, .4, .4);
	// Right
    glBegin(GL_QUADS);
		glNormal3d(0, 0, 1);
		glTexCoord2f(0, 0);
		glVertex3dv(c);
		glTexCoord2f(1, 0);
		glVertex3dv(g);
		glTexCoord2f(1, 1);
		glVertex3dv(f);
		glTexCoord2f(0, 1);
		glVertex3dv(b);
    glEnd();

    glDisable(GL_TEXTURE_2D);
}


void CParticleSystem::drawGlass()
{
	float bx = (maxX+1)/2;
	float bz = (maxZ+1)/2;

	float height = maxY + .25;

	//Radius of particles = .2
	//glColor4d(.68, .81, .88, .5);
	float a[] = {-bx, (float)cy, -bz};
	float b[] = {bx-pRadius, (float)cy, -bz};
	float c[] = {bx-pRadius, (float)cy, bz};
	float d[] = {-bx, (float)cy, bz};

	GLfloat black[] = {.0, .0, .0, 1.0};
	GLfloat white[] = {1.0, 1.0, 1.0, 1.0};
	glMaterialfv(GL_FRONT, GL_AMBIENT, black);
	glMaterialfv(GL_FRONT, GL_DIFFUSE, black);
	glLightfv(GL_LIGHT0, GL_SPECULAR, white);
	glLightfv(GL_LIGHT1, GL_SPECULAR, white);

	glDisable(GL_LIGHTING);
	glEnable (GL_BLEND);
	glBlendFunc (GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

	
	float s = .15f; //size of the sides
	float aL[] = {a[0]+s, a[1], a[2]};
	float aR[] = {a[0], a[1], a[2]+s};
	float bL[] = {b[0], b[1], b[2]+s};
	float bR[] = {b[0]-s, b[1], b[2]};
	float cL[] = {c[0]-s, c[1], c[2]};
	float cR[] = {c[0], c[1], c[2]-s};
	float dL[] = {d[0], d[1], d[2]-s};
	float dR[] = {d[0]+s, d[1], d[2]};

	glMaterialfv(GL_FRONT, GL_DIFFUSE, white);

	glPushMatrix();
	glTranslatef(bx, 0, bz);

	//Draw edges of box first
	glColor4f(.8f, .8f, .8f, 0.5f);
	drawCorner(a, aL, aR, height, 0);
	drawCorner(b, bL, bR, height, 0);
	drawCorner(c, cL, cR, height, 0);
	drawCorner(d, dL, dR, height, 0);
	
	drawCorner(a, aL, aR, maxZ+1., 1);
	drawCorner(b, bL, bR, maxZ+1., 1);
	drawCorner(a, aL, aR, maxX+.5, 2);
	drawCorner(d, dL, dR, maxX+.5, 2);
	
	float aH[] = {a[0], a[1] + maxY + .25, a[2]};
	float bH[] = {b[0], b[1] + maxY + .25, b[2]};
	float dH[] = {d[0], d[1] + maxY + .25, d[2]};
	float aHL[] = {a[0]+s, a[1]+maxY + .25, a[2]};
	float aHR[] = {a[0], a[1]+maxY + .25, a[2]+s};
	float bHL[] = {b[0], b[1]+maxY + .25, b[2]+s};
	float bHR[] = {b[0]-s, b[1]+maxY + .25, b[2]};
	float dHL[] = {d[0], d[1]+maxY + .25, d[2]-s};
	float dHR[] = {d[0]+s, d[1]+maxY + .25, d[2]};
	drawCorner(aH, aHL, aHR, maxZ+1., 1);
	drawCorner(bH, bHL, bHR, maxZ+1., 1);
	drawCorner(aH, aHL, aHR, maxX+.5, 2);
	drawCorner(dH, dHL, dHR, maxX+.5, 2);

	//Now the glass panes
	glColor4f(.8f, .8f, .8f, 0.1f);
	glBegin(GL_QUADS);
		glVertex3f(d[0],d[1]+height,d[2]);
		glVertex3fv(d);
		glVertex3fv(a);
		glVertex3f(a[0],a[1]+height,a[2]);
	glEnd();
	glBegin(GL_QUADS);
		glVertex3f(a[0],a[1]+height,a[2]);
		glVertex3fv(a);
		glVertex3fv(b);
		glVertex3f(b[0],b[1]+height,b[2]);
	glEnd();
	glBegin(GL_QUADS);
		glVertex3fv(b);
		glVertex3f(b[0],b[1]+height,b[2]);
		glVertex3f(c[0],c[1]+height,c[2]);
		glVertex3fv(c);
	glEnd();
	glBegin(GL_QUADS);
		glVertex3fv(c);
		glVertex3f(c[0],c[1]+height,c[2]);
		glVertex3f(d[0],d[1]+height,d[2]);
		glVertex3fv(d);
	glEnd();
	glBegin(GL_QUADS);
		glVertex3f(b[0],b[1]+height,b[2]);
		glVertex3f(a[0],a[1]+height,a[2]);
		glVertex3f(d[0],d[1]+height,d[2]);
		glVertex3f(c[0],c[1]+height,c[2]);
	glEnd();

	glDisable(GL_BLEND);
	glPopMatrix();
	
}

void drawCorner(float *cor, float *L, float *R, float ht, int dir)
{
	float corH[] = {cor[0], cor[1] + ht, cor[2]};
	float LH[] = {L[0], L[1] + ht, L[2]};
	float RH[] = {R[0], R[1] + ht, R[2]};

	if (dir == 1) {
		corH[1] -= ht, corH[2] += ht;
		LH[1] -= ht, LH[2] += ht;
		RH[1] -= ht, RH[2] += ht;
	}
	else if (dir == 2) {
		corH[1] -= ht, corH[0] += ht;
		LH[1] -= ht, LH[0] += ht;
		RH[1] -= ht, RH[0] += ht;
	}


	glBegin(GL_QUAD_STRIP);
		glVertex3fv(cor);
		glVertex3fv(corH);
		glVertex3fv(R);
		glVertex3fv(RH);
		glVertex3fv(L);
		glVertex3fv(LH);
		glVertex3fv(cor);
		glVertex3fv(corH);
	glEnd();
}


void CParticleSystem::drawSphereAt(Vector3df loc, GLfloat c[3], bool alpha) 
{
	GLfloat white[] = {1.f, 1.f, 1.f, 1.0f};

	float x= 0.f, y= 0.f, z= 0.f;
	x = loc.x();
	y = loc.y();
	z = loc.z();

	glPushMatrix();
	glTranslatef(x, y, z);
	if (alpha) {
		glMaterialfv(GL_FRONT, GL_DIFFUSE, white);
		glColor4f(c[0], c[1], c[2], 0.2f);	
	}
	else {
		glMaterialfv(GL_FRONT, GL_DIFFUSE, white);
		glColor4f(c[0], c[1], c[2], 1.f);	
	}

	//Approximation of simple line tetrahedron -- alternate to drawing spheres
	/*float sz = pRadius;
	glBegin(GL_LINE_LOOP);
		glVertex3f(0.f,0.f,0.f);//a
		glVertex3f(sz,0.f,0.f);//b
		glVertex3f(0.f,sz,0.f);//c
		glVertex3f(0.f,0.f,0.f);//a
		glVertex3f(0.f,0.f,sz);//d
		glVertex3f(0.f,sz,0.f);//c
		//glVertex3f(sz,0.f,0.f);//b //Don't really need all the lines for appx
		//glVertex3f(0.f,0.f,sz);//d
	glEnd();*/

	//Quadrics are part of glu library that make parameterized shapes
    GLUquadricObj *qobj;
    qobj = gluNewQuadric();
    gluQuadricDrawStyle(qobj, GLU_FILL); 

	gluSphere(qobj, pRadius, 20, 10);
	glPopMatrix();
}




//--------------------------------Creating the setting-------------------------------
//-----------------------------------------------------------------------------------

const GLfloat  AMBIENT[] = {0.1f, 0.1f, 0.1f, 1.0f};
const GLfloat  LIGHT0LOC[]    = {1.f, 0.5f, 1.f, 0.f};
const GLfloat  LIGHT0DCOLOR[] = {1.f, 1.f, 1.f, 1.f};
const GLfloat  LIGHT0SCOLOR[] = {1.f, 1.f, 1.f, 1.f};
const GLfloat  LIGHT1LOC[]    = {-2.f, 0.4f, 1.f, 0.f};
const GLfloat  LIGHT1DCOLOR[] = {1.f, 1.f, 1.f, 1.f};
const GLfloat  LIGHT1SCOLOR[] = {1.f, 1.f, 1.f, 1.f};
const GLfloat  LIGHT2LOC[]    = {0.f, -2.f, -2.f, 0.f};
const GLfloat  LIGHT2DCOLOR[] = {1.f, 1.f, 1.f, 1.f};
const GLfloat  LIGHT2SCOLOR[] = {1.f, 1.f, 1.f, 1.f};

void CParticleSystem::SetLighting()
{
    GLfloat black[] = {0.f, 0.f, 0.f, 1.0f};

    glLightModelfv(GL_LIGHT_MODEL_AMBIENT, black);

    glEnable(GL_LIGHT0);
    glLightfv(GL_LIGHT0, GL_POSITION, LIGHT0LOC);
    glLightfv(GL_LIGHT0, GL_DIFFUSE, LIGHT0DCOLOR);
    glLightfv(GL_LIGHT0, GL_SPECULAR, LIGHT0SCOLOR);
    glLightfv(GL_LIGHT0, GL_AMBIENT, AMBIENT);

    glEnable(GL_LIGHT1);
    glLightfv(GL_LIGHT1, GL_POSITION, LIGHT1LOC);	
    glLightfv(GL_LIGHT1, GL_DIFFUSE, LIGHT1DCOLOR);
    glLightfv(GL_LIGHT1, GL_SPECULAR, LIGHT1SCOLOR);

	/*glEnable(GL_LIGHT2);
	glLightfv(GL_LIGHT2, GL_POSITION, LIGHT2LOC);
    glLightfv(GL_LIGHT2, GL_DIFFUSE, LIGHT2DCOLOR);
    glLightfv(GL_LIGHT2, GL_SPECULAR, LIGHT2SCOLOR);*/
}


//Helper function, here for now
float randFloat(float low, float high)
{
       //Source: www.daniweb.com/software-development/cpp/threads/98545
       return (float)(rand() / (static_cast<float>(RAND_MAX) + 1.0)) * (high - low) + low;
}
